/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"

#include "volFields.H"
#include "fvcGrad.H"

// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //

template<class Type>
void Foam::functionObjects::tkeFluxBudget::processField
(
    const word& fieldName,
    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tvalue
)
{
    typedef GeometricField<Type, fvPatchField, volMesh> FieldType;

    const word scopedName = word("tke_") + fieldName;

    if (obr_.foundObject<FieldType>(scopedName))
    {
        FieldType& fld =
            const_cast<FieldType&>(obr_.lookupObject<FieldType>(scopedName));
        fld == tvalue();
    }
    else if (obr_.found(scopedName))
    {
        WarningInFunction
            << "Cannot store turbulence field " << scopedName
            << " since an object with that name already exists"
            << nl << endl;
    }
    else
    {
      obr_.store
        (
            new FieldType
            (
                IOobject
                (
                    scopedName,
                    obr_.time().timeName(),
                    obr_,
                    IOobject::READ_IF_PRESENT,
                    IOobject::NO_WRITE
                ),
                tvalue
            )
        );
    }
}


// returns the convection term
Foam::tmp<Foam::volScalarField>
Foam::functionObjects::tkeFluxBudget::Ck() const
{
  const volVectorField& U = obr_.lookupObject<volVectorField>("U");

  // trace of Reynolds stress tensor
  const volSymmTensorField& UP2M = obr_.lookupObject<volSymmTensorField>("UPrime2Mean");

  // inner product of uprime
  const volScalarField ipUPrime( tr(UP2M) );
  
  // gradient of tIPUPrime
  const volVectorField gradIPUPrime(fvc::grad(ipUPrime));

  // inner product of U and tgradIPUPrime
  const volScalarField ipUUPrime( U & gradIPUPrime );

  return tmp<volScalarField>
  (
      new volScalarField
      (
          IOobject
	  (
	      "Ck",
	      ipUUPrime.mesh().time().timeName(),
	      ipUUPrime.mesh()
	   ),
	  (-0.5)*ipUUPrime,
	  ipUUPrime.boundaryField().types()
       )
   );

}


// returns the production term
Foam::tmp<Foam::volScalarField>
Foam::functionObjects::tkeFluxBudget::Pk() const
{
  const volVectorField& U = obr_.lookupObject<volVectorField>("U");
  const volSymmTensorField& UP2M = obr_.lookupObject<volSymmTensorField>("UPrime2Mean");

  // gradient of velocity
  const volTensorField gradU( fvc::grad(U) );

  // inner product of UP2M and gradU
  const volScalarField ipUP2MGradU( UP2M && gradU );

  return tmp<volScalarField>
  (
      new volScalarField
      (
          IOobject
	  (
	      "Pk",
	      ipUP2MGradU.mesh().time().timeName(),
	      ipUP2MGradU.mesh()
	   ),
	  (-1.0)*ipUP2MGradU,
	  ipUP2MGradU.boundaryField().types()
       )
   );
}


// returns the turbulence transport term
Foam::tmp<Foam::volScalarField>
Foam::functionObjects::tkeFluxBudget::Tk() const
{
  const volVectorField& U = obr_.lookupObject<volVectorField>("U");
  const volVectorField& UMean = obr_.lookupObject<volVectorField>("UMean");
  
  // fluctutation field
  const volVectorField UPrime = U - UMean;

  // inner product of UPrime and UPrime
  const volScalarField ipUPrime( UPrime & UPrime );

  // outer product of fluctuation field
  const volVectorField opUPrime( ipUPrime * UPrime );

  // divergence of opUPrime
  const volScalarField divOpUPrime( fvc::div(opUPrime) );


//add turbulent fluctuation heat flux


//add turbulent fluctuation heat flux 2021.3.29
  return tmp<volScalarField>
  (
      new volScalarField
      (
          IOobject
	  (
	      "Tk",
	      divOpUPrime.mesh().time().timeName(),
	      divOpUPrime.mesh()
	   ),
	  (-0.5)*divOpUPrime,
	  divOpUPrime.boundaryField().types()
       )
   );

}


// returns the viscous diffusion term
Foam::tmp<Foam::volScalarField>
Foam::functionObjects::tkeFluxBudget::Dk() const
{
  const volSymmTensorField& UP2M = obr_.lookupObject<volSymmTensorField>("UPrime2Mean");

  // trace of UP2M
  const volScalarField trUP2M( tr(UP2M) );

  // laplacian of trUP2M
  const volScalarField lapUP2M( fvc::laplacian(trUP2M) );

  return tmp<volScalarField>
  (
      new volScalarField
      (
          IOobject
	  (
	      "Dk",
	      lapUP2M.mesh().time().timeName(),
	      lapUP2M.mesh()
	   ),
	  nu_tkeFlux*(0.5)*lapUP2M,
	  lapUP2M.boundaryField().types()
       )
   );
}


// returns the viscous dissipation term
Foam::tmp<Foam::volScalarField>
Foam::functionObjects::tkeFluxBudget::Epik() const
{
  const volVectorField& U = obr_.lookupObject<volVectorField>("U");
  const volVectorField& UMean = obr_.lookupObject<volVectorField>("UMean");
 
  // fluctuation field
  const volVectorField UPrime = U - UMean;

  // gradient of the fluctuation field
  const volTensorField gradUPrime( fvc::grad(UPrime) );

  // double inner product of gradUPrime
  const volScalarField ipGradUPrime( gradUPrime && gradUPrime );
  
  return tmp<volScalarField>
  (
      new volScalarField
      (
          IOobject
	  (
	      "Epik",
	      ipGradUPrime.mesh().time().timeName(),
	      ipGradUPrime.mesh()
	   ),
	  (-1.0)*nu_tkeFlux*ipGradUPrime,
	  ipGradUPrime.boundaryField().types()
       )
   );
}


// returns the velocity-pressure gradient term
Foam::tmp<Foam::volScalarField>
Foam::functionObjects::tkeFluxBudget::Pik() const
{
  const volVectorField& U = obr_.lookupObject<volVectorField>("U");
  const volVectorField& UMean = obr_.lookupObject<volVectorField>("UMean");
 
  const volScalarField& P = obr_.lookupObject<volScalarField>("p");
  const volScalarField& PMean = obr_.lookupObject<volScalarField>("pMean");

  // fluctuation field
  const volVectorField UPrime = U - UMean;
  const volScalarField PPrime = P - PMean;

  // gradient of the PPrime field
  const volVectorField gradPPrime( fvc::grad(PPrime) );

  // double inner product of UPrime and gradPPrime
  const volScalarField ipUPrimeGradPPrime( UPrime & gradPPrime );
  
  return tmp<volScalarField>
  (
      new volScalarField
      (
          IOobject
	  (
	      "Pik",
	      ipUPrimeGradPPrime.mesh().time().timeName(),
	      ipUPrimeGradPPrime.mesh()
	   ),
	  (-1.0)*ipUPrimeGradPPrime,
	  ipUPrimeGradPPrime.boundaryField().types()
       )
   );
}

// ************************************************//

//1//
Foam::tmp<Foam::volVectorField>
Foam::functionObjects::tkeFluxBudget::Hf_0_001() const
{
  const volVectorField& U = obr_.lookupObject<volVectorField>("U");
  const volVectorField& UMean = obr_.lookupObject<volVectorField>("UMean");
 
  const volScalarField& T_0_001 = obr_.lookupObject<volScalarField>("T_0_001"); 
  const volScalarField& T_0_001Mean = obr_.lookupObject<volScalarField>("T_0_001Mean");

  // fluctuation field
  const volVectorField UPrime = U - UMean;
  const volScalarField T_0_001Prime = T_0_001 - T_0_001Mean;

  const volVectorField heatfluxPrime_0_001(UPrime * T_0_001Prime);
  return tmp<volVectorField>
  (
      new volVectorField
      (
          IOobject
	  (
	      "Hf_0_001",
	      heatfluxPrime_0_001.mesh().time().timeName(),
	      heatfluxPrime_0_001.mesh()
	   ),
	  (-1.0)*heatfluxPrime_0_001,
	  heatfluxPrime_0_001.boundaryField().types()
       )
   );
}

//2//
Foam::tmp<Foam::volVectorField>
Foam::functionObjects::tkeFluxBudget::Hf_0_005() const
{
  const volVectorField& U = obr_.lookupObject<volVectorField>("U");
  const volVectorField& UMean = obr_.lookupObject<volVectorField>("UMean");
 
  const volScalarField& T_0_005 = obr_.lookupObject<volScalarField>("T_0_005"); 
  const volScalarField& T_0_005Mean = obr_.lookupObject<volScalarField>("T_0_005Mean");

  // fluctuation field
  const volVectorField UPrime = U - UMean;
  const volScalarField T_0_005Prime = T_0_005 - T_0_005Mean;

  const volVectorField heatfluxPrime_0_005(UPrime * T_0_005Prime);
  return tmp<volVectorField>
  (
      new volVectorField
      (
          IOobject
	  (
	      "Hf_0_005",
	      heatfluxPrime_0_005.mesh().time().timeName(),
	      heatfluxPrime_0_005.mesh()
	   ),
	  (-1.0)*heatfluxPrime_0_005,
	  heatfluxPrime_0_005.boundaryField().types()
       )
   );
}

//3//
Foam::tmp<Foam::volVectorField>
Foam::functionObjects::tkeFluxBudget::Hf_0_025() const
{
  const volVectorField& U = obr_.lookupObject<volVectorField>("U");
  const volVectorField& UMean = obr_.lookupObject<volVectorField>("UMean");
 
  const volScalarField& T_0_025 = obr_.lookupObject<volScalarField>("T_0_025"); 
  const volScalarField& T_0_025Mean = obr_.lookupObject<volScalarField>("T_0_025Mean");

  // fluctuation field
  const volVectorField UPrime = U - UMean;
  const volScalarField T_0_025Prime = T_0_025 - T_0_025Mean;

  const volVectorField heatfluxPrime_0_025(UPrime * T_0_025Prime);
  return tmp<volVectorField>
  (
      new volVectorField
      (
          IOobject
	  (
	      "Hf_0_025",
	      heatfluxPrime_0_025.mesh().time().timeName(),
	      heatfluxPrime_0_025.mesh()
	   ),
	  (-1.0)*heatfluxPrime_0_025,
	  heatfluxPrime_0_025.boundaryField().types()
       )
   );
}
//4//
Foam::tmp<Foam::volVectorField>
Foam::functionObjects::tkeFluxBudget::Hf_0_05() const
{
  const volVectorField& U = obr_.lookupObject<volVectorField>("U");
  const volVectorField& UMean = obr_.lookupObject<volVectorField>("UMean");
 
  const volScalarField& T_0_05 = obr_.lookupObject<volScalarField>("T_0_05"); 
  const volScalarField& T_0_05Mean = obr_.lookupObject<volScalarField>("T_0_05Mean");

  // fluctuation field
  const volVectorField UPrime = U - UMean;
  const volScalarField T_0_05Prime = T_0_05 - T_0_05Mean;

  const volVectorField heatfluxPrime_0_05(UPrime * T_0_05Prime);
  return tmp<volVectorField>
  (
      new volVectorField
      (
          IOobject
	  (
	      "Hf_0_05",
	      heatfluxPrime_0_05.mesh().time().timeName(),
	      heatfluxPrime_0_05.mesh()
	   ),
	  (-1.0)*heatfluxPrime_0_05,
	  heatfluxPrime_0_05.boundaryField().types()
       )
   );
}

//5//
Foam::tmp<Foam::volVectorField>
Foam::functionObjects::tkeFluxBudget::Hf_0_125() const
{
  const volVectorField& U = obr_.lookupObject<volVectorField>("U");
  const volVectorField& UMean = obr_.lookupObject<volVectorField>("UMean");
 
  const volScalarField& T_0_125 = obr_.lookupObject<volScalarField>("T_0_125"); 
  const volScalarField& T_0_125Mean = obr_.lookupObject<volScalarField>("T_0_125Mean");

  // fluctuation field
  const volVectorField UPrime = U - UMean;
  const volScalarField T_0_125Prime = T_0_125 - T_0_125Mean;

  const volVectorField heatfluxPrime_0_125(UPrime * T_0_125Prime);
  return tmp<volVectorField>
  (
      new volVectorField
      (
          IOobject
	  (
	      "Hf_0_125",
	      heatfluxPrime_0_125.mesh().time().timeName(),
	      heatfluxPrime_0_125.mesh()
	   ),
	  (-1.0)*heatfluxPrime_0_125,
	  heatfluxPrime_0_125.boundaryField().types()
       )
   );
}

//6//
Foam::tmp<Foam::volVectorField>
Foam::functionObjects::tkeFluxBudget::Hf_0_25() const
{
  const volVectorField& U = obr_.lookupObject<volVectorField>("U");
  const volVectorField& UMean = obr_.lookupObject<volVectorField>("UMean");
 
  const volScalarField& T_0_25 = obr_.lookupObject<volScalarField>("T_0_25"); 
  const volScalarField& T_0_25Mean = obr_.lookupObject<volScalarField>("T_0_25Mean");

  // fluctuation field
  const volVectorField UPrime = U - UMean;
  const volScalarField T_0_25Prime = T_0_25 - T_0_25Mean;

  const volVectorField heatfluxPrime_0_25(UPrime * T_0_25Prime);
  return tmp<volVectorField>
  (
      new volVectorField
      (
          IOobject
	  (
	      "Hf_0_25",
	      heatfluxPrime_0_25.mesh().time().timeName(),
	      heatfluxPrime_0_25.mesh()
	   ),
	  (-1.0)*heatfluxPrime_0_25,
	  heatfluxPrime_0_25.boundaryField().types()
       )
   );
}

//7//
Foam::tmp<Foam::volVectorField>
Foam::functionObjects::tkeFluxBudget::Hf_0_71() const
{
  const volVectorField& U = obr_.lookupObject<volVectorField>("U");
  const volVectorField& UMean = obr_.lookupObject<volVectorField>("UMean");
 
  const volScalarField& T_0_71 = obr_.lookupObject<volScalarField>("T_0_71"); 
  const volScalarField& T_0_71Mean = obr_.lookupObject<volScalarField>("T_0_71Mean");

  // fluctuation field
  const volVectorField UPrime = U - UMean;
  const volScalarField T_0_71Prime = T_0_71 - T_0_71Mean;

  const volVectorField heatfluxPrime_0_71(UPrime * T_0_71Prime);
  return tmp<volVectorField>
  (
      new volVectorField
      (
          IOobject
	  (
	      "Hf_0_71",
	      heatfluxPrime_0_71.mesh().time().timeName(),
	      heatfluxPrime_0_71.mesh()
	   ),
	  (-1.0)*heatfluxPrime_0_71,
	  heatfluxPrime_0_71.boundaryField().types()
       )
   );
}

//8//
Foam::tmp<Foam::volVectorField>
Foam::functionObjects::tkeFluxBudget::Hf_1() const
{
  const volVectorField& U = obr_.lookupObject<volVectorField>("U");
  const volVectorField& UMean = obr_.lookupObject<volVectorField>("UMean");
 
  const volScalarField& T_1 = obr_.lookupObject<volScalarField>("T_1"); 
  const volScalarField& T_1Mean = obr_.lookupObject<volScalarField>("T_1Mean");

  // fluctuation field
  const volVectorField UPrime = U - UMean;
  const volScalarField T_1Prime = T_1 - T_1Mean;

  const volVectorField heatfluxPrime_1(UPrime * T_1Prime);
  return tmp<volVectorField>
  (
      new volVectorField
      (
          IOobject
	  (
	      "Hf_1",
	      heatfluxPrime_1.mesh().time().timeName(),
	      heatfluxPrime_1.mesh()
	   ),
	  (-1.0)*heatfluxPrime_1,
	  heatfluxPrime_1.boundaryField().types()
       )
   );
}



// ************************************************************************* //
